﻿# -*- coding: utf-8 -*-
from copy import deepcopy
from random import randint

class Matrix:
    """
    A naïve (and native) Python implementation of matrices with basic operations.
    by Guillaume HUYSMANS, October 2013 (LGPL license)
    
    Benchmarks (made with linearEquations using this library):
    * xxx
    * yyy
    """
    
    def __init__(self, data):
        """data may be an existing Matrix object to copy OR an array"""
        if isinstance(data, Matrix):
            self.arr = deepcopy(data.arr)
            self.rows = data.rows
            self.cols = data.cols
        else:
            #It must be an array...
            c = 0
            for r in data:
                if len(r)>c:
                    c = len(r)
                elif len(r)<c:
                    raise ValueError("The matrix is uncomplete!")
            self.cols = c
            self.rows = len(data)
            if self.rows==0:
                raise(ValueError("A matrix can't have no rows!"))
            if self.cols==0:
                raise(ValueEror("A matrix can't have an empty row"))
            self.arr = deepcopy(data)
        
    def __str__(self):
        """Returns a string representation of the matrix (no alignment)"""
        s = ''
        for r in self.arr:
            for c in r:
                s += str(c)+','
            s = s[:-1]+'\n'
        return s[:-1]

    def isSquare(self):
        """Returns a boolean indicating whether the matrix is square"""
        return self.cols==self.rows

    def subMatrix(self, r, c):
        """Returns a new matrix without row r and column c"""
        res = deepcopy(self.arr)
        if r<self.rows:
            res.pop(r) #remove the row
        if c<self.cols:
            for r in res:
                r.pop(c) #remove the column
        return Matrix(res)

    def trans(self):
        """Computes the transposed matrix"""
        if self.rows==1:
            return Matrix([[v] for v in self.arr[0]])
        else:
            return Matrix(zip(*self.arr))

    def cofact(self):
        """Computes the cofactors matrix"""
        res = []
        row = [0]*self.cols
        for i, r in enumerate(self.arr):
            for j, v in enumerate(r):
                row[j] = ((-1)**(i+j)) * self.subMatrix(i,j).det()
            res.append(deepcopy(row))
        return Matrix(res)

    def det(self):
        """
        Computes the matrix's determinant (it must be square).
        It uses an hard-coded formula for 2x2 matrices and
        the cofactors method for larger ones.
        """
        if self.rows != self.cols:
            return #undefined!
        elif self.rows==1:
            return self.arr[0][0] #extreme case...
        elif self.rows==2:
            #needed by cofact (faster than calling det on a 1x1 matrix)
            return self.arr[0][0]*self.arr[1][1]-self.arr[1][0]*self.arr[0][1]
        else:
            res = 0
            for i, v in enumerate(self.arr[0]):
                res += ((-1)**i) * self.subMatrix(0,i).det() * v
            return res

    def getCol(self, c):
        """Returns a list containing values from the c-th column"""
        res = [0] * self.rows
        for i in range(self.rows):
            res[i] = self.arr[i][c]
        return res

    def setCol(self, c, data):
        """
        Replaces the c-th column by the values given in data
        data may be a n*1 matrix or a list
        """
        if isinstance(data, Matrix):
            if data.cols != 1:
                raise(TypeError("data must be a 1*n matrix!"))
            data = [v[0] for v in data.arr]
        for i, v in enumerate(data):
            self.arr[i][c] = data[i]

    def __eq__(self, other):
        """Compares two matrices (without epsilon)"""
        return self.arr == other.arr

    def __add__(self, other):
        if isinstance(other, Matrix) == False:
            raise(TypeError("Can't add something else than a matrix to a matrix!"))
        if self.rows!=other.rows or self.cols!=other.cols:
            raise(TypeError('Both matrix must be the same size!'))
        res = []
        for i, r in enumerate(self.arr):
            res.append([sum(t) for t in zip(r, other.arr[i])])
        return Matrix(res)

    def __neg__(self):
        return self*(-1)

    def __sub__(self, other):
        return self*(-1)+other

    def __mul__(self, other):
        if isinstance(other, Matrix):
            if self.cols!=other.rows:
                return
            res = [[0]*other.cols for x in range(self.rows)]
            #res = Matrix.null(self.rows, other.cols) #FIXME why doesn't it work?
            #For each cell in the result array,
            for i in range(self.rows):
                for j in range(other.cols):
                    #Compute the dot product between a row from self
                    #and a column from other and store it
                    for a in range(self.cols):
                        res[i][j] += self.arr[i][a]*other.arr[a][j]
            return Matrix(res)
        else:
            return Matrix([[float(other)*v for v in r] for r in self.arr])

    def __div__(self, other):
        return self*(1/float(other))

    def __rdiv__(self, other):
        #FIXME fails with matA/matA!!
        #other must be a number: if it were a matrix,
        #this method wouldn't get called.
        inv = self.cofact().trans()/self.det()
        #2) multiply the result with other
        return inv*other

    def __pow__(self, exp):
        """Raises the matrix to the exp-th power (exp must be an integer)"""
        if isinstance(exp, Matrix):
            raise(TypeError("The exponent can't be a matrix, too!"))
        exp = int(exp)
        #Get a positive exponent (and invert the matrix if needed)
        if exp<0:
            res = 1/self
            exp *= -1
        else:
            res = Matrix(self) #copy
        #Multiply as much as necessary
        while exp>1:
            res *= self
            exp -= 1
        return res

    def __iadd__(self, other):
        return self+other

    def __isub__(self, other):
        return self-other

    def __imul__(self, other):
        return self*other

    def __ipow__(self, exp):
        return self**exp

    def __radd__(self, other):
        return self+other #will always cause an error (other must be a number)
        #this function is actually useless (just here to customize the exception)

    def __rmul__(self, other):
        #other must be a number: if it were a matrix,
        #this method wouldn't get called.
        return self*other

    def iden(n):
        """Returns an identity matrix (n*n)"""
        res = []
        row = [0] * n
        for i in range(n):
            if i>0:
                row[i-1] = 0
            row[i] = 1
            res.append(deepcopy(row))
        return Matrix(res)
    iden = staticmethod(iden)

    def null(r, c):
        """Returns an r*c matrix filled with zeroes"""
        #FIXME unusable in __mul__ --> why?
        return Matrix([[0]*r for x in range(c)])
    null = staticmethod(null)

    def random(r, c, x=-100, y=100):
        """Returns an r*c matrix filled with random values between x and y"""
        return Matrix([[randint(x,y) for n in range(c)] for o in range(r)])
    random = staticmethod(random)



#Testing stub
if __name__ == "__main__":
    matA = Matrix([[2,3,1],[4,5,2],[3,2,-1]])
    matB = Matrix([[5,4,2],[6,7,6],[2,3,10]])
    id3 = Matrix.iden(3)
    id3.setCol(0, [9, 9, 9])
    print "setC1\n",id3,"\n"
    print "getC1\n",id3.getCol(0),"\n"
    print "A\n",matA,"\n"
    print "A^(-1)\n",1/matA,"\n"
    print "A * 1/A = identity3\n",matA*(1/matA),"\n"
    print "1/A * A = identity3\n",(1/matA)*matA,"\n"
    print "AB\n",matA*matB,"\n"
    print "BA (different)\n",matB*matA,"\n"
